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Abstract 

C^ A simple algorithm and a computational program to numerically compute the electric field gradient and 

^ ^ the concomitant quadrupolar nuclear splitting is developed for an arbitrary ionic crystal. The calculations 

Ch are performed using a point charge model. The program provides three different ways for the data input: 

, ^ by Bravais lattices, by lattice parameters, or by introducing any spatial structure. The program calculates 

^.^ the components of the electric field gradient, the asymmetry parameter and the quadrupolar splitting for a 

~: given number of nearest neighbors with respect to the nuclear charge as origin. In addition, the program 
allows the use of different Sternheimer antishielding factors. 

* r^ Keywords: Electric field gradient, quadrupolar splitting, Mossbauer spectroscopy, algorithm and 

C/3 numerical computation, asymmetry parameter, crystallographic lattices 
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^ 1. Introduction 

■^ 

C^ The electrostatic energy W due to the interaction of a nuclear charge distribution p{r) and the electrostatic 

C potential V{r) generated by its electric environment is given by 

'§ W= f p{i)V{i)d\ (1) 

^S. J vol 

where d^r is the volume element and r = {xi^x2,x^) are spatial coordinates. The integral is calculated 
over the nucleus volume. A suitable way to evaluate it is to make a multipole expansion of the electrostatic 
potential V{r) around the center of charge of the nucleus as origin, assuming that V{x) is a slowly varying 



> 

in 



f-^ function over the nuclear region. Expanding in a Taylor serie around the nucleus center of charge, one 

\f^ obtains CEl 
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^—s The relevant terms in this expansion are the first and third terms, due to the fact that the second one 

,-H is zercp] because, when multiplied by the nuclear charge, it represents the interaction of the nuclear dipole 

^""l moment (w^hich is zero) with the external electric field, E. The next non-zero terms are several orders 

J> of magnitude smaller than the third one-^J so, in a good approximation, the interaction energy can be 

expressed as 

3 



X 

3 vv = - 

C3 2 



^ ^ = J E ^1^^^^ (3) 



j,fe= 



* Corresponding author: jorge_hdz@ciencias.unam.mx 
** Principal corresponding author: rgomez@servidor.unam.mx, tel: (+52) (55) 5622-4849 
^ As well as all the odd terms in the expansion. 



Preprint submitted to Elsevier July 4, 2012 



where Vjk are the electric field gradient (EFG) tensor components, and Qij are the quadrupolar nuclear 
moment components; both are second rank tensors. Choosing a principal axis system for the EFG tensor, 
the interaction energy can be expressed as the sum of two terms 
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The first term, called isomer shift, represents the effect due to the nucleus size. The second one corresponds 
to the so called quadrupolar nuclear splitting AQ, so the interaction hamiltonian between the nuclear 
quadrupolar moment Q and the EFG tensor WE, with respect to an arbitrary axes system with origin in 
the nuclear charge centroid, is given by 

'k = --eQ(E)WE (5) 
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where ® denotes tensorial product. Considering the z axis along the largest component of the EFG (Vzz = 
eg) and the Laplace equation, the hamiltonian |p| is transformed through the Wigner-Eckart theorem 1^ into 

where P, Ix, ly and I^ are the nuclear spin magnitudes, and 

V = 7^-^^ (7) 

is the so called asymmetry parameter, which indicates how much the electric potentml departs from spherical 
symmetry. An analytical solution of Ipl can only be obtained for the / = ■^/2 case ' ^' ' . By far, the most used 
isotope in Mossbauer spectroscopy is Fe, for which the useful transition is I — ^/2 ^ I — 1/2 and, in what 
follows, we will restrict to this case. The analytical solution is: 



^ e qQ 2 _ ^2] L _ rf (8) 
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That is, the nuclear / — 2/2 energy level is split into two levels (±2/2 and ±1/2) and the ground level 
/ = 1/2 stays degenerated. This gives rise to two absorption lines in the Mossbauer spectrum separated by 
an energy 

AQ = ^/^ (9) 

which is called the quadrupolar nuclear splitting. 

1.1. Electric Gradient Tensor 

In rectangular coordinates, the EFG of a set of n point charges is: 

^^,^.=Eg^-( '^'";f""" ) (10) 

where qk and r^. = (2:1^ , 2:2^ , 2^3^ ) are the respective charge and position of the /c*'* ion. The electric interaction 
of the nucleus with its surroundings has two different origins: the charge density of the electrons of the 
nucleus under study, and the ligands of the crystal lattice™!. 

For ionic crystals, the main contribution to the EFG comes from those ions directly coordinated to the 
nucleus. Considering that the interatomic distances in a crystal are pretty much larger than the displacements 
due to the ion vibrations, a good approximation to the EFG can be done using equation ( [10| with a point 
charge model. The contribution to the EFG due to the electronic distribution requires a complex calculation 



that involves not only the knowledge of the electronic wave functions of the atom, but also shielding and 
antishielding effects and polarization of the electronic distribution due to other charges near by. This type 
of calculations retreat from the scope of this work, so their effect will be taken into account trough two 
parameters R and 700, known as the shielding and antishielding Stemheimer factors'^ ^^ , to obtain 

q = qugi'^ - loo) + qvaii'^ ~ R) (11) 

The first term in equation lITTJ corresponds to the ligand contribution, while the second one constitutes 
the valence contribution. However, in ionic crystals the last one can be neglected. The following program 
is developed for this case. 

2. Structure of the computational program 

The program was focused as an useful tool in a Mossbauer spectroscopy laboratory, so it computes the 
components of the EFG tensor and the quadrupolar splitting for a ^^Fe nucleus by default. However, it is 
able to work with any nucleus, just by introducing its respective quadrupolar moment value. 

In order to compute the EFG in a great number of crystalline lattices, three different input data modes 
were developed, allowing for a wide range of applications. Those modes are briefly described here: 

• N arbitrary ions: In this section of the program, the coordinates and valences of each ligand constituting 
the crystalline array are inputted manually. The spatial distribution of the ions can be totally arbitrary. 
The algorithm can handle a number of ions as large as necessary, being this number, of course, finite. 

• Bravais lattices: Here, an election of one of the fourteen possible Bravais Lattices in three dimensions 
is made, just by introducing the parameter(s) that define such lattice. The program allows to select 
the place in the lattice in w^hich the EFG will be computed. 

• Lattice parameters: When available, one can input the values of the lattice parameters, so the program 
identifies the respective lattice, reconstructing it in order to carry out the computations. 

The program was developed in a structured computational language, so there is a main module calling 
different functions and subroutines. 

2.1. Functions 

There are three functions defined in the program. 

2.1.1. R 

The function R calculates the euclidean rectangular distance between the «*'' ligand coordinates and the 
nucleus under study taken as the origin. 

2.1.2. V 

This function calculates the Vx^Xi component of the EFG tensor in principal axes for the i*'' ion with 
equation 1 10 1, for each value of k, and where Qk is the valence charge of the ligand, qug. 



2.1.3. DQ 

This function computes the value of the quadrupolar splitting as a function of Vzz and the asymmetry 
parameter, through equations (|9l and | [TT| , leaving the result in terms of the (1 — 700) factor, without 
considering the valence contribution in the total charge. 



2.2. Main module 

Here, the physical constants to be used in the program are defined. It establishes the value of the 
quadrupolar moment Q to be usecFl and the ionic configuration to be worked out, which are: 

2.2.2. N arbitrary ions 

• Introduce the number N of ligands to be considered in the computation. 

• Introduce the three coordinates (in angstroms) and the valence of each ion. 

• The distance to the origin is computed for each ion, via the R function. 

• The components of the EFG are calculated, adding^ to each one the contribution of each ion through 
function V. 

• The largest component of the EFG is assigned to \Vzz\, and \Vyy\ > \Vxx\- 

• The asymmetry parameter is computed with equation |l7||. 

• The value of the quadrupolar splitting is computed with the function DQ. 

• The results of the EFG components, the asymmetry parameter and the quadrupolar splitting are 
shown on screen and saved in a file. 

2.2.2. Bravais Lattices 

• Choose one of the seven possible groups in three dimensions. 

• Select the lattice to be taken into accounrland introduce the parameter(s) that define it. Then choose 
the number of nearest neighbors to be deemed, the valence of each layer of neighbors, and the position 
in the structure in w^hich the EFG is to be computed (in the center or the vertex of the structure). 

• With this information, the program reckons the coordinates of the ligands in the lattice, through the 
algorithm presented in the next section. Once the coordinates are determined, the distance of each 
ligand to the origin is computed. In order to identify and count the layers of nearest neighbors, 
the information of all the generated neighbors is ordered and displayed in growing distances to the 
origirrl The components of the EFG are calculated for the chosen layers of neighbors. 

• The largest component of the EFG assigned to \Vzz\, and \Vyy\ > \ Vxx \ ■ 

• The asymmetry parameter is computed with equation |l7||. 

• The value of the quadrupolar splitting is reckoned with the function DQ. 

• The results of the EFG components, the asymmetry parameter and the quadrupolar splitting are 
shown on screen and saved in a file. 

2.2.3. Lattice parameters 

• Introduce the six lattice parameters, (a, h, c) and (a, (5, 7), in angstroms and degrees respectively. 

• The program identifies the lattice that corresponds with the lattice parameters introduced, and with 
the information of the lattice, the program proceeds in the same form than in the previous section, 

Bravais lattices. 



^The program uses the value Q = 16b of the quadrupolar nuclear moment for the ^^Fe, recently reported by Martinez-Pinedo et 
alEH. 

^Due to the superposition principle of the electric potential. 

*If possible, choose between the simple, the body centered, the two faces centered or the face centered correspondent structure. 

*The ordering process is carried out with a bubble ordering algorithm, which does not compromises the efficiency of the program 
because the lists of neighbors to be ordered are not generally too large in standard calculations in solid state. 



3. Algorithm 

In what follows, the main algorithm used by the sections Bravais lattices and lattice parameters to find the 
points in the lattice where the ligands are to be considered for the computations, is described: 

a) Select the number of nearest neighbors to be considered. 

b) • If all the neighbors have the same valence, introduce it. 

• If not, introduce the valence of each layer of neighbors. 

c) Choose to compute the EFG in the center or in the vertex of the structure. 

d) With the six lattice parameters (a, b, c) y {a, (3, 7), the rectangular components of the crystallographic 



axes are calculated through the next transformation equations, obtained in the appendix Appendix A 



ax ^ a (12) 

bx = 6 cos 7 
by — 6 sin 7 
Cx = c cos /3 
Cy = c(cos a CSC 7 — cos /? cot 7) 
Cz = c sin /3-\/l — (cos a esc /3 esc 7 — cot /3 cot 7)^ 
e) If the studied nucleus is centered in the body, the coordinates of the ligands are calculated as follows: 
• If the structure is simple (SC, ST, SO, SM, triclinic, trigonal or hexagonalrl the coordinates of 



the i*'* ion are computed through equations 1 13 ', w^here the numbers ni, 712 and n^ are whole 
numbers in the interval [— m, mYj 

Xi = {ni + -)ax + (n2 + -)bx + (ng + -)cx (13) 

Zi = (na + -)cz 

• If the structure is body centered (BCC, BCT or BCOp, the coordinates of the i*'* ion are computed 
through equations | [14| , the same way as in simple structures, but excluding the point (0, 0, 0). 

Xi ^ niax + n2bx + n-iCx (14) 

Hi = n2by + n^Cy 
Zi = nzCz 

^These are common abbreviations in crystallography. SC: Simple Cubic; ST: Simple Tetragonal; SO: Simple Orthorhombic; SM: 
Simple Monoclinic; BCC: Body Centered Cubic; BCT: Body Centered Tetragonal; BCO: Body Centered Orthorhombic; FCC: Face 
Centered Cubic; FCO: Face Centered Orthorhombic; 2FCO: Two Face Centered Orthorhombic; 2FCM: Two Face Centered Monoclinic. 

^The number m (chosen in the step a)) depends on the number of nearest neighbors to be considered in the computation. For 
example, m = 2 is enough to find the third nearest neighbors. 



If the structure is face centered (FCC or FCOp, the coordinates of the i*^ ion are computed 
through equations | [T5) , 
planes, where the numb 
ni y 71.2 can not be zero. 



through equations u5\, u6\ and 1 17 1 for the ligands in the faces parallel to the crystallographic 
planes, where the numbers ni, rt2 y «3 are w^hole numbers in the interval [—m, mPand such that 



Xi = -nia^ + -n2bx + n^Cx (15) 
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= »^3C2 



2/i = ;:;"-2&a + n'iCy 



-- -^nia^ + -n2Cx + n^bx (16) 

y* = 7,'^2Cy + n^by 
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2-2^ 



a^^j = x'^iCi: + 7;'^2bx + n^ax (17) 

1 , 1 
Vi = ijn2by + -riiCy 

1 

Zi = l^'^lCz 



If the structure is two face centered (2FCO or 2FCM^, the coordinates of the i*-^ ion are computed 
through equations 1 18 i, w^here the numbers ni, n2 y n^ are whole numbers in the interval [~m, m^ 



and such that n^ can not be zero. 

Xt = niax + n2bx + ^n^Cx (18) 



X + n2bx + -n^c 


n2by + -nsCy 


1 



f) If the studied nucleus is centered in the vertex of the structure, the coordinates of the ligands are 
computed as follows: 

• If the structure is simple (SC, ST, SO, SM, triclinic, trigonal or hexagonal)^, the coordinates of the 
I*'' ion are computed through equations | [l4| , excluding the point (0, 0, 0). 

• If the structure is body centered (BCC, BCT or BCOp, the coordinates of the i*'* ion are computed 
through equations | [T9| , the same way as in simple structures, excluding the point (0, 0, 0). 

Xi = -niUx + -n2bx + -n^Cx (19) 

1 , 1 
Vi = ^n2by + -WiCy 

1 

Zi = o"^'^-' 

• If the structure is face centered (FCC or FCOp, the coordinates of the i*'' ion are computed 



through equations | [15| , | [16| and 1 17 1 for the ligands in the faces parallel to the crystallographic 
planes, where the numbers ni, n2 and 713 are whole numbers in the interval [— m, mPand such 
that ni y ri2 can not be zero. 



If the structure is two face centered (2FCO or 2FCMpi the coordinates of the i*'' ion are computed 
through equations 1 20 i, where the numbers ni, n2 y n^ are whole numbers in the interval [—to, to-P 



and such that ni and n2 can not be zero. 



^nia^ + -n2b^ + n^c^ (20) 



1 
2 



-n2by + n-iCy 



g) The distance to the origin is computed for each ion, via the R function, 
h) The ligands are ordered in growing distances to the origirP. 
i) The number of ions in each layer of the nearest neighbors is counted^ 
j) The valences introduced in the step b) are assigned to each layer of the nearest neighbors computed 

previously, 
k) The components of the EFG are computed, adding^ to each one the contribution of each ion through 
function V. 

It is important to point out that in our calculations the contribution of the Sternheimer factors to the EFG 
has not been taken into account, so the calculated values will be normally smaller than the experimental 
values. However, the purpose of these calculations is to endow a guide to discriminate between different 
site environments of the iron nucleus through the relative magnitudes of the calculated EFG ^ . 

4. Conclusions 

The algorithm and the program developed here are useful as a high applicability tool in both experimental 
spectroscopy and in any theoretical research in solid state and crystallography, requiring this kind of 
computations. 

The program presented in this work is extremely versatile and friendly with the final user, and only 
requires the structural information of the system under study. In spite of the fact that the EFG tensor and 
its quadrupolar splitting computations are based in a simple point charge model, disregarding the valence 
contribution, it is a useful tool to discriminate the different structures present in complex ionic systems. 

The program was written in Fortran 77, to assure high compatibility across different platforms, but 
the structure of the program and the algorithm are equally useful if the program is written in any other 
structured computational language, like C, Phyton, Pascal, etc., without compromising its accuracy and 
stability. 

Of course, the program can be improved including the effect of the electronic density in the EFG 
tensor, but that matter may be explored in a future work. However, the computation of the shielding and 
antishielding factors is complicated, so their effect has to be taken into account with empirical adjustments, 
without changing significantly the relative magnitudes of the quadrupolar splittings. 

Appendix A. Transformation from crystallographic to rectangular coordinates 



Consider the crystallographic axes a, b and c, and the rectangular ones x, y and z, as shown in figure 
The crystallographic axes can be expressed in the cartesian base, as: 



A.l 



*" Indeed, the program can be used only to count the number of neighbors, their coordinates and distances to the origin for a wide 
range of lattices. 




Figure A.l: Geometrical scheme of the coordinate 
transformation 



a = aCx 

b ~ b cos 763; + b sin 761, 

c = f cos TiCx + f sin TiCy + hcz 



(A.l) 
(A.2) 
(A.3) 



Solving the system of equations that arise from the seven triangles shown in figure A.l and considering 
the respective constrictions, a general expression of the crystallographic axes in the rectangular base in 
terms of the lattice parameters can be written as 

a — acx (A.4) 

b — bcos^ex + bsin^fey (A.5) 

c = c cos /3e^ + c(cos a esc 7 — cos /3 cot 7)6^ (A.6) 

+ c sin /3 yjl — (cos a esc /3 esc 7 — cot (3 cot 7)^6^ 
for any non-orthonormal set of crystallographic axes. 
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